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1 . SUMMARY 


This is the final summary report by Resource Management Associates, 
Inc., of the first year's work under Contract No. NASW-4011 to the 
National Aeronautics and Space Administration. The work under this 
initial phase of the contract relates to the preliminary development of 
a global, 3-D magnetohydrodynamic computational model to quantitatively 
describe the detailed continuum field and plasma interaction process of 
the solar wind with cometary and planetary bodies throughout the solar 
system. The work extends a highly-successful , observationally-verif ied 
computational model previously developed by the author, and is 
appropriate for the global determination of supersonic, ' super-Alfvenic 
solar wind flows past planetary obstacles. This report provides a 
concise description of the problems studied, a summary of all the 
important research results, and copies of the publications. 
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2. PROGRAM OVERVIEW 


The overall objective of the present program relates to the 
development of a global, 3-D magnetohydrodynamic computational model to 
quantitatively describe the detailed continuum field and plasma 
interaction process of the solar wind with cometary and planetary 
bodies. The work involves detailed extensions of a highly-successful , 
observational ly-verified solar wind/terrestrial planet interaction model 
previously developed by the principal investigator. That model, for 
reasons of computational efficiency, was preliminarily developed based 
upon a gas dynamic plus convected field approximation to the full MHD 
equations. In this program that model is being extended to the general 
3-D MHD level and will additionally be generalized to include mass 
addition effects due to such phenomena as photoionization and charge 
exchange. The technical work of this study involes new computational 
development and numerical modeling of a variety of detailed and advanced 
aspects of the 3-D solar wind interaction problem. The ultimate goal is 
the development of a comprehensive, user-oriented, solar 
wind/cometary-planetary obstacle interaction computational model that 
would embody a hierarchy of three different accuracy levels and 
associated computational algorithms. For reasons that were described in 
detail in the original proposal, this hierachy of different accuracy 
levels is: 


Level 


Model 


1 Axisymmetric gas dynamic + convected B field 

2 3-D gas dynamic + 3-D convected B field 

3 3-D MHD 


A major objective of the first phase of the present study is to 
develop the Level 2 3-D gas dynamic plus 3-D convected B field model. 
For reasons also described in detail in the original proposal, all 3 
levels of the model will contain two separate flow field solvers that 
are coupled together to determine the steady state flow field. This 
feature of our model, which . is unique among all currently existing 
global interaction models, allows very high resolution of the entire 
flow field, particularly in the tail region. This is illustrated in 
Figure 1, where a typical flow field grid density of the combined Level 
1 axisymmetric flow solvers is shown. The region from the subsolar 
point to the terminator plane is solved via a time marching solver (NOSE 
code) to the steady state; whereas the region downstream of the 
terminator is solved via a spatially marching solver (TAIL code) which 
advances the solution downstream as far as required. As can .be seen, 
the downstream grid density remains at very high resolution. Such 
detail of the downstream flow field cannot be achieved by any of the 
other currently existing interaction models which employ a single 
unsteady flow solver to determine the entire nose plus tail flow 
regions. The major technical tasks involved in the development of the 
Level 2 model are the development of the 3-0 versions of the NOSE and 
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TAIL flow field codes (N0SE3D, TAIL3D). In addition, an entirely new 
magnetic field solver (MAG3D) to determine the convected B field in a 
general 3-D gas dynamic flow field is also required. This requirement 
is necessary since the decomposition method of Alksne and Webster 
(Ref. 1) employed in the Level 1 axisymmetric model is valid only for 
axisymmetric flow fields and cannot be applied to general 3-D flows. 

In addition to the development of the 3-D Level 2 model, several 
other technical tasks were planned for this first phase. These involve 
several extensions of the Level 1 model, and include preliminary 
development of a mass loading capability -to account for such phenomena 
as photoionization and charge exchange, a complete diagnostic and 
analysis graphics package for detailed presentation of plasma and field 
properties throughout the magnetosheath, and a self-consistent procedure 
to determine the correct rather than approximate magnetoionopause shape 
and which locally satisfies the appropriate balance between 
magnetosheath and magnetospheric or ionospheric pressures across the 
magnetoionopause boundary. Finally, in addition to these model 
developments, a variety of collaborative investigations and applications 
involving the model were undertaken by us together with a number of 
space scientists. All of these developments are discussed in the 
following section. 
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3. PHASE I ACCOMPLISHMENTS 


The major technical accomplishment of the Phase I effort involved 
the development of the sophisticated 3-D Level 2 gas dynamic flow field 
and convected B field solvers. To accomplish this, it was necessary to 
develop and integrate three separate codes. Two of these codes are for 
determination of the flow field, and correspond to the 3-D versions of 
the NOSE and TAIL codes (N0SE3D, TAIL3D) shown in Figure 1. The third 
code (MAG3D) determines the convected B field in a general 3-D gas 
dynamic flow field. We provide a brief . discussion of the development of 
these codes below. 

After an extensive review of existing methods, selection of the 
basic solver to be employed in the N0SE3D code was made. The 
methodology is based upon the AIR3D code developed at NASA/Ames Research 
Center (Ref. 2). The solution algorithm employed is the Beam and 
Warming implicit approximate factorization algorithm used in delta form 
described in Ref. 3. The original basic algorithm is first-order 
accurate in time, noniterative, second-order accurate in spatial 
derivatives with central difference operators. This methodology is 
capable of solving both the inviscid gas dynamic Euler equations as well 
as, if required, a certain form of the viscous gas dynamic equations 
known as the thin-layer Navier-Stokes equations. These equations 
include certain additional terms over and above the Euler equation level 
which represent viscous effects in the body-normal direction. This 
capability is important to our study for future planned applications in 
which we wish to study certain viscous effects near the magnetoionopause 
boundary. The computational mesh employed in this method is one which 
incorporates a fitted shock outer boundary which exactly satisfies the 
unsteady, 3-D Rankine-Hugoniot relations and moves with the shock during 
the convergence process. The inner boundary is implemented as a fitted 
stationary surface representing the magnetoionopause surface. A mesh 
clustering capability for high density clustering in the body-normal 
direction has been implemented specifically to enable accuracy to be 
maintained for viscous calculations. 

The gas dynamic solver employed for the Level 2 solution of the 3-D 
tail region of the flow field is based on the 3-D gas dynamic marching 
code reported in Ref. 4. The solution algorithm is a shock capturing, 
fully conservative form that employs second-order noncentered spatial 
derivatives to solve the steady, inviscid, gas dynamic Euler equations. 
As with. the N0SE3D code, the mesh to be used with the TAIL3D code will 
employ a fitted outer boundary for the bow shock, at which location the 
steady 3-D Rankine-Hugoniot relations are exactly satisfied. The inner 
boundary of the mesh will also be fitted but to an impenetrable surface 
representing the magnetoionopause boundary. 

The method chosen to determine the convected B field for a general 
3-D gas dynamic flow field involves the 3-D generalization of. the 
kinematic procedure described by Spreiter in Ref. 5. In our 
implementation of the method, we have formulated a particularly 
efficient procedure. We initiate a streamline computation at the point 
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in the magnetosheath where the B field solution is desired. The 
streamline computation is then carried upstream through the bow shock 
into the uniform oncoming flow. At a selected point on that streamline 
in the uniform oncoming flow, a vector of incremental length centered 
about that point is formed in the direction of the oncoming IMF. Two 
streamlines that bracket the original streamline and which past through 
the ends of the incremental g vector are then integrated downstream, at 
equal time steps with the original streamline, until they reach the 
original point in the magnetosheath where the B field is desired. The 
direction and magnitude of the B field are then given directly by the 
stretched and rotated incremental B vector. This procedure is repeated 
at each point in the magnetosheath where the 3-D magnetic field is 
required. 

Thus far in this phase, we have completed the initial development 
of all three of the above codes. The N0SE3D code, which requires by far 
the most computational resources since it integrates the unsteady 3-D 
Euler equations forward in time to a steady state, has now been 
developed and verified on a number of benchmark cases. We have 
similarly completed development and verification of the TAIL3D code. 
Finally, we have coupled the N0SE3D and TAIL3D codes via a joining data 
plane at the terminator location, and have initially verified the 
coupling procedure by performing several benchmark cases in which TAIL3D 
solutions were carried sufficiently far downstream from the terminator 
to demonstrate that a stable and accurate coupled solution procedure has 
been achieved. Similarly, we have developed the MAG3D code, and have 
initially verified it via comparisons with the Level 1 B field procedure 
for. the special case of axisymmetric flow. 

In addition to the Level 2 3-D model development, we have developed 
in this phase a number of extensions and enhancements in the Level 1 
model. These include a complete diagnostic and analysis graphics 
package that provides all of the detailed plasma and field properties 
throughout the interactive magnetosheath region. The package includes: 
maps of the plasma streamlines, velocity, temperature, and density 
contours; magnetic field contours and field line patterns; spacecraft 
trajectory diagrams in ecliptic and polar plane view; time histories of 
plasma density, temperature, velocity, magnitude and three velocity 
components; magnetic field . magnitude and three magnetic field 
components; maps of entire computational grid; plots displaying the 
variation along arbitrary segments of the computational mesh of velocity 
magnitude and three velocity components, density, pressure, temperature, 
and magnetic field magnitude and three magnetic field components. The 
graphics package was intentionally based on CALCOMP routines for easy 
portability to other computer facilities. All of the plots included in 
the package have now been verified, and the package has been developed 
for use. on a high-density dot matrix printer, a feature never before 
available with the model. The output from the graphics package is 
essentially report quality, and the capability provides for the first 
time the means for instantaneous viewing of the model results at almost 
no cost. Figures 2-4 provide some examples of an abbreviated subset of 
the complete grahical output of the analysis package, while Figures 5-7 
provide corresponding examples of the diagnostic package. 
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With regard to the preliminary development of a mass loading 
capability in the Level 1 model, at this time we have completed 
derivation of the appropriate source model representation for mass 
loading due to photoionoization, and have also derived the appropriate 
modifications to the partial differential equations describing the 
conservation of mass, momentum, energy, and magnetic flux when such 
source terms due to this phenomena are present. Final development, 
individual implementation into both the NOSE and TAIL codes, and 
verification of the capability, together with case studies at Venus, 
will be completed in Phase II. 

We have also proceeded to develop a self-consistent iterative 
procedure in the Level 1 model to improve the satisfaction of the 
pressure balance boundary condition at the magnetoionopause boundary. 
The present Level 1 model employs a specified magnetoionopause shape 
that is previously determined using a Newtonian surface pressure 
variation from the stagnation point to points along the surface up to 
the terminator. This variation is used to approximate the pressure 
distribution on the magnetosheath side of the magnetoionopause. The 
improved procedure we are developing will iteratively alter that initial 
shape so as to ultimately satisfy the true pressure balance that locally 
exists between the exact magnetosheath plasma and magnetospheric/iono- 
spheric pressure at each point along the magnetoionopause boundary. 
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4. COLLABORATIVE INVESTIGATIONS 


Finally, in addition to all of the model developments described 
above, we have also carried out a variety of collaborative 
investigations involving use of the model. These studies were 
undertaken by us together with a number of space scientists. We 
consider these investigatons to be one of the most important aspects of 
the model development being undertaken here, as they effectively 
guarantee the widespread utility of the model through interdisciplinary 
studies involving many different topics of current space science 
interest. 

A list of the publications that have resulted from those particular 
studies and which are now complete is as follows: 


Publications from Completed Collaborative Studies 

Magnetopause Merging Site Asymmetries, N. U. Crooker, J. G. Luhmann, 
J. R. Spreiter and S. S. Stahara. Jour. Geophys. Res., Vol . 90, No. Al, 
January 1985, pp. 341-346. 

Magnetic Field Draping Against the Dayside Magnetopause, N. U. Crooker,. 
J. G. Luhmann, C. T. Russell, E. J. Smith, J. R. Spreiter and S. S. Sta- 
hara. Jour. Geophys. Res., Vol. 90, No. A4, April 1985, pp. 3505-3510. 


Solar Wind Flow About the Outer Planets: Gasdynamic Modeling of the 
Jupiter and Saturn Bow Shocks, J. A. Slavin, E. J. Smith, J. R. Spreiter 
and S. S. Stahara. Jour. Geophys. Res., Vol. 90, No. A7, July 1985, 


pp. 6275-6286. 

Evidence for . Mass-Loading of the 
C. T. Russell, J. R. Spreiter and S. 
5, No. 4, 1985, pp. 307-311. 

A Gas Dynamic Magnetosheath Field 
Fields: Application to the Solar 
J. G. Luhmann, R. J. Warniers, C. 
S. S. Stahara. Jour. Geophys. Res., 
pp. 3001-3010. 

Magnetohydrodynamic and Gasdynamic 
J. R. Spreiter and S. S. Stahara. 
Collisionless Shock Waves in the 
Research, 1985. 1985. 


Venus Magnetosheath, J. G. Luhmann, 
S. Stahara. Adv. Space Res., Vol. 


Model for Unsteady Interplanetary 
Wind Interaction with Venus, 
T. Russell, J. R. Spreiter and 
Vol. 91, No. A3, March 1985, 


Theories for Planetary Bow Waves, 
Geophysical Monograph 35, 
Heliosphere: Reviews of Current 


In addition to these completed studies, we are also currently 
engaged in other collaborative studies involving use of the model. 
These include: (1) studies involving the quantitative determination of 
3-D effects of the solar wind interaction with the asymmetric 
magnetopauses of Jupiter and Saturn - with Dr. James Slavin and 
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Dr. Edward Smith of JPL; (2) studies of the predicted and observed 
ionopause shapes at Venus when the exact local plasma pressure in the 
magnetosheath is balanced with the observed local ionospheric pressure - 
with Dr. William Knudsen of Knudsen Geophysical Research. 
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Figure 1. Illustration of typical flow-field grid density for 
Level 1 model; = 3.0, y = 5/3, H/R 0 = 0.03 
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Figure 2. Sample flow field output from Level 1 model graphical output package: 
density contours for n 3.0, y M 5/3 flow past an ionopause shape 
with H/R 0 = 0.03 
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Figure 4. Sample time history output from Level 1 model graphical 
output package: temperature vs. time for a trajectory 

passing through the magnetosheath region and ionopause 
of the flow field shown in Figure 1 
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Figure 5. Sample plot from Level 1 model extended diagnostic graphic output package: 
computational grid for M ro = 3.0, Y= 5/3 flow past an ionopause shape with 


VELOCITY 



Figure 6 . Sample plasma property output from Level 1 model extended diagnostic graphic output 
package: velocity vector plots along constant radial mech locations at one mesh 

spacing away from the obstacle surface (1=2) and along bow shock (1=19) 



UNIT MAGNETIC FIELD : PERPENQICULAR COMPONENT 



Figure 7. Sample magnetic field output from Level 1 model extended diagnostic graphic 
output package: vector plots of unit magnetic field perpendicular component 

along various radial mcch linos normal to the obstacle for M - 3.0, Y s 5 

flow past an ionopause shape with H/R q = 0.03 
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